Install/load necessary R packages
Set date range
min_year = 2017
max_year = 2021
More information about peatland and upland water elevation data collected at the Marcell Experimental Forest can be found on the Environmental Data Initiative Repository Portal.
Read in the most up-to-date data from the following Box file
path:
External-MEF_DATA/Hydro/WaterTableElevation/Bogwells/L1
External-MEF_DATA/EDI/peatland_water_table_daily/edi.562.2/data_objects/Clean the data into a tidy format
#set file path to read in data from Box
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/Hydro/WaterTableElevation/Bogwells/L1"
#read in the raw data
bogwell <- read_csv(here(filepath, "L1_Peatland_well_daily_history_2021forMia.csv"))
#clean the data
##note: we do not want any E values for the flag column
bogwell_clean <- bogwell %>%
clean_names() %>%
mutate(year = format(as.POSIXct(date, format = '%Y-%m-%d',
tz = "GMT"),
format = '%Y'),
year = as.numeric(year)) %>%
subset(year >= min_year & year <= max_year)
#save clean data CSV to use in future RMD files
write.csv(x = bogwell_clean,
file = file.path(here("intermediate_data", "bogwell_clean.csv")),
row.names = FALSE)
Plot all sites together
#plot
p_bogwell <- ggplot(data = bogwell_clean) +
geom_line(aes(x = date, y = wte, color = peatland)) +
theme_classic() +
labs(x = "Time",
y = "Water Table Elevation (m above sea level)",
title = paste0("Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")"),
color = "Site"
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5,
size = 7))
#static plot
p_bogwell
bogwell_clean_s2 <- bogwell_clean %>%
subset(peatland == "S1")
#plot
p_bogwell_s1 <- ggplot(data = bogwell_clean_s2) +
geom_line(aes(x = date,
y = wte,
group = 1,
text = paste0("Date: ", date, "\n",
"Water Table Elevation (m): ", wte))) +
theme_classic() +
labs(x = "Time",
y = "Water Table Elevation (m above sea level)",
title = paste0("S1 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")")
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5,
size = 7))
#interactive plot
ggplotly(p_bogwell_s1,
tooltip = "text")
bogwell_clean_s2 <- bogwell_clean %>%
subset(peatland == "S2")
#plot
p_bogwell_s2 <- ggplot(data = bogwell_clean_s2) +
geom_line(aes(x = date,
y = wte,
group = 1,
text = paste0("Date: ", date, "\n",
"Water Table Elevation (m): ", wte))) +
theme_classic() +
labs(x = "Time",
y = "Water Table Elevation (m above sea level)",
title = paste0("S2 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")")
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5,
size = 7))
#interactive plot
ggplotly(p_bogwell_s2,
tooltip = "text")
bogwell_clean_s2 %>%
#remove rows that have NA wte values
filter(!is.na(wte)) %>%
summarise(mean_wte_m = round(mean(wte), digits = 3),
max_wte_m = max(wte),
min_wte_m = min(wte),
sd_wte_m = round(sd(wte), digits = 3)
) %>%
#t() %>%
kable(caption = "S2 Stripchart Bogwell Statistics",
col.names = c("Mean WTE (m)",
"Max WTE (m)",
"Min WTE (m)",
"SD")) %>%
kable_material(c("striped", "hover"))
| Mean WTE (m) | Max WTE (m) | Min WTE (m) | SD |
|---|---|---|---|
| 421.88 | 422.15 | 421.48 | 0.095 |
bogwell_clean_s3 <- bogwell_clean %>%
subset(peatland == "S3")
#plot
p_bogwell_s3 <- ggplot(data = bogwell_clean_s3) +
geom_line(aes(x = date,
y = wte,
group = 1,
text = paste0("Date: ", date, "\n",
"Water Table Elevation (m): ", wte))) +
theme_classic() +
labs(x = "Time",
y = "Water Table Elevation (m above sea level)",
title = paste0("S3 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")")
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5,
size = 7))
#interactive plot
ggplotly(p_bogwell_s3,
tooltip = "text")
bogwell_clean_s4 <- bogwell_clean %>%
subset(peatland == "S4")
#plot
p_bogwell_s4 <- ggplot(data = bogwell_clean_s4) +
geom_line(aes(x = date,
y = wte,
group = 1,
text = paste0("Date: ", date, "\n",
"Water Table Elevation (m): ", wte))) +
theme_classic() +
labs(x = "Time",
y = "Water Table Elevation (m above sea level)",
title = paste0("S4 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")")
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5,
size = 7))
#interactive plot
ggplotly(p_bogwell_s4,
tooltip = "text")
bogwell_clean_s5 <- bogwell_clean %>%
subset(peatland == "S5")
#plot
p_bogwell_s5 <- ggplot(data = bogwell_clean_s5) +
geom_line(aes(x = date,
y = wte,
group = 1,
text = paste0("Date: ", date, "\n",
"Water Table Elevation (m): ", wte))) +
theme_classic() +
labs(x = "Time",
y = "Water Table Elevation (m above sea level)",
title = paste0("S5 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")")
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5,
size = 7))
#interactive plot
ggplotly(p_bogwell_s5,
tooltip = "text")
bogwell_clean_s6 <- bogwell_clean %>%
subset(peatland == "S6")
#plot
p_bogwell_s6 <- ggplot(data = bogwell_clean_s6) +
geom_line(aes(x = date,
y = wte,
group = 1,
text = paste0("Date: ", date, "\n",
"Water Table Elevation (m): ", wte))) +
theme_classic() +
labs(x = "Time",
y = "Water Table Elevation (m above sea level)",
title = paste0("S6 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")")
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5,
size = 7))
#interactive plot
ggplotly(p_bogwell_s6,
tooltip = "text")
Function for loading Campbell DataLogger .dat files:
readCDL = function(file){
# read data file starting on 5th line
dat <- read.csv(file, sep=",",header=FALSE,skip=4,na.strings="NAN",stringsAsFactors=F)
# Read in just the header line (l2)
# unlist the line, and remove quotes
h <- readLines(file, n=2)[2]
n <- as.factor(unlist(strsplit(h, ",")) )
n2 <- gsub('"', "", n)
# assign column names to dataframe
colnames(dat) = n2
return(dat)
}
Constants to convert between feet and meters:
MetersPerFoot = 0.3048
FeetPerMeter = 3.28048
min_year = 2020
max_year = 2021
Read in the 2022 data from the following Box file path:
External-MEF_DATA/DataDump/RealTimeData
Clean the 2022 data into a tidy format
#set file path to read in data from Box
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/RealTimeData"
#read in the raw data
##note: this data is only 2022
s1_se <- readCDL(file = here(filepath, "S1-EM3_Table1.dat"))
#clean the data
s1_se_clean <- s1_se %>%
clean_names() %>%
#select relevant columns
select(timestamp,
wt_elev) %>%
#rename columns
rename(datetime = timestamp,
wt_elev_ft = wt_elev) %>%
#reformat columns into appropriate classes
mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S',
tz = "GMT"),
wt_elev_ft = as.numeric(wt_elev_ft),
wt_elev_m = wt_elev_ft * MetersPerFoot)
2019 data
#set file path to read in data from Box
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/annual_appended_logger_files"
year = "2019"
file = "S1-EM3_Table1.dat"
#read in 2019 data
s1_2019 <- readCDL(file = here(filepath, year, file))
#clean 2019 data
s1_2019_clean <- s1_2019 %>%
clean_names() %>%
#select relevant columns
select(timestamp,
wt_elev) %>%
#rename columns
rename(datetime = timestamp,
wt_elev_ft = wt_elev) %>%
#reformat columns into appropriate classes
mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S',
tz = "GMT"),
wt_elev_ft = as.numeric(wt_elev_ft),
wt_elev_m = wt_elev_ft * MetersPerFoot)
2020 data
#set file path to read in data from Box
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/annual_appended_logger_files"
year = "2020"
file = "S1-EM3_Table1.dat"
#read in 2020 data
#note: this data is 2020 - 2022
s1_2020 <- readCDL(file = here(filepath, year, file))
#clean 2020 data
s1_2020_clean <- s1_2020 %>%
clean_names() %>%
#select relevant columns
select(timestamp,
wt_elev) %>%
#rename columns
rename(datetime = timestamp,
wt_elev_ft = wt_elev) %>%
#reformat columns into appropriate classes
mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S',
tz = "GMT"),
wt_elev_ft = as.numeric(wt_elev_ft),
wt_elev_m = wt_elev_ft * MetersPerFoot)
Combine S1 data
s1_all <- rbind(s1_se_clean, s1_2019_clean, s1_2020_clean) %>%
#remove NA values
filter(!is.na(wt_elev_m)) %>%
#extract year
mutate(year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S',
tz = "GMT"),
format = '%Y'),
year = as.numeric(year))
#set min and max year
min_year = min(s1_all$year)
max_year = max(s1_all$year)
p_s1 <- ggplot(data = s1_all) +
geom_line(aes(x = datetime,
y = wt_elev_m,
group = 1,
text = paste0("Date: ", datetime, "\n",
"Water Table Elevation (m): ", round(wt_elev_m, digits = 3)))) +
theme_classic() +
labs(x = "Time",
y = "Water Table Elevation (m above sea level)",
title = paste0("S1 Shaft Encoder Water Table Elevation (", min_year, "-", max_year, ")")
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5,
size = 7))
#interactive plot
ggplotly(p_s1,
tooltip = "text")
MaxWTElev column#set file path to read in data from Box
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/RealTimeData"
#read in the raw data
s2_se <- readCDL(file = here(filepath, "S2_bogwell_S2BW_Daily.dat"))
s2_se_clean <- s2_se %>%
clean_names() %>%
#select relevant columns
select(timestamp,
max_wt_elev) %>%
#rename columns
rename(datetime = timestamp,
max_wt_elev_ft = max_wt_elev) %>%
#reformat columns into appropriate classes
mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S',
tz = "GMT"),
max_wt_elev_ft = as.numeric(max_wt_elev_ft),
max_wt_elev_m = max_wt_elev_ft * MetersPerFoot)
#plot
p_s2 <- ggplot(data = s2_se_clean) +
geom_line(aes(x = datetime,
y = max_wt_elev_m,
group = 1,
text = paste0("Date: ", datetime, "\n",
"Water Table Elevation (m): ", max_wt_elev_m))) +
theme_classic() +
labs(x = "Time",
y = "Water Table Elevation (m)",
title = paste0("S2 Bogwell Water Table Elevation (", min_year, "-", max_year, ")")
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
#interactive plot
ggplotly(p_s2,
tooltip = "text")
s2_se_clean %>%
#rename columns
rename(wte = max_wt_elev_m) %>%
#remove rows that have NA wte values
filter(!is.na(wte)) %>%
summarise(mean_wte_m = round(mean(wte), digits = 3),
max_wte_m = max(wte),
min_wte_m = min(wte),
sd_wte_m = round(sd(wte), digits = 3)
) %>%
#t() %>%
kable(caption = "S2 Shaft Encoder Bogwell Statistics",
col.names = c("Mean WTE (m)",
"Max WTE (m)",
"Min WTE (m)",
"SD")) %>%
kable_material(c("striped", "hover"))
| Mean WTE (m) | Max WTE (m) | Min WTE (m) | SD |
|---|---|---|---|
| 421.844 | 422.0916 | 421.4854 | 0.115 |
#format stripchart data
bog <- bogwell_clean_s2 %>%
#remove unnecessary columns
select(-peatland, -flag) %>%
#create a column to note what type of data this is
mutate(type = "stripchart")
#format shaft encoder data
se <- s2_se_clean %>%
#extract date and year
mutate(date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S',
tz = "GMT"),
format = '%Y-%m-%d'),
date = as.POSIXct(date, tz = "GMT"),
year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S',
tz = "GMT"),
format = '%Y'),
year = as.numeric(year),
#note data type
type = "shaft_encoder") %>%
#remove unnecessary columns
select(-datetime, -max_wt_elev_ft) %>%
#rename columns to join the data
rename(wte = max_wt_elev_m)
#join the data
s2_join <- full_join(x = bog,
y = se,
by = c("wte", "date", "year", "type"))
#subset the join to the time period where the two data sources overlap
s2_join_sub <- s2_join %>%
pivot_wider(names_from = "type",
values_from = "wte") %>%
#remove rows with NA values
filter(!is.na(shaft_encoder),
!is.na(stripchart))
p <- ggplot(data = s2_join_sub) +
geom_point(aes(x = stripchart,
y = shaft_encoder,
alpha = 0.1,
group = 1,
text = paste0("Stripchart: ", stripchart, "\n",
"Shaft Encoder: ", shaft_encoder))) +
#plot straight line to compare point spread
geom_abline(intercept = 0,
slope = 1) +
theme_classic() +
labs(x = "Manual Stripchart",
y = "Shaft Encoder",
title = paste0("S2 1:1 WTE Comparison (", min_year, "-", max_year, ")")
) +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5,
size = 7),
legend.position = "none")
#interactive plot
ggplotly(p = p,
tooltip = "text")
Note that this file is being knit as index.html
into the climate_bogwell
repository in order to update the GitHub pages
website, where the plots can be viewed online.